Chapter 8 HMSC analysis

8.1 Load data

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

8.2 Variance partitioning

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()

variable	mean	sd
Random: location	37.900015	25.317903
logseqdepth	56.110626	25.796874
sex	4.937460	5.612719
origin	1.051899	1.282563

# Basal tree
varpart_tree <- genome_tree

#Varpart table
varpart_table <- varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
varpart_tree <- varpart_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
   scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()

# Add variance stacked barplot
vertical_tree <-  varpart_tree +
       scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
        geom_fruit(
             data=varpart_table,
             geom=geom_bar,
             mapping = aes(x=value, y=genome, fill=variable, group=variable),
             pwidth = 2,
             offset = 0.05,
             width= 1,
             orientation="y",
             stat="identity")+
      labs(fill="Variable")

vertical_tree

8.3 Posterior estimates

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    #select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
    column_to_rownames(var="genome")

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)


colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap

postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

postestimates_tree +
        vexpand(.25, 1) # expand top 

8.3.1 Positively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()

genome	phylum	class	order	family	species	value
bin_CaboVerde.50	Actinomycetota	Actinomycetia	Actinomycetales	Actinomycetaceae	NA	0.953
bin_Aruba.25	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium pseudocatenulatum	0.995
bin_Aruba.2	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium adolescentis	0.994
bin_Aruba.6	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium longum	0.985
bin_Denmark.14	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium gallinarum	0.950
bin_Aruba.47	Actinomycetota	Actinomycetia	Mycobacteriales	Mycobacteriaceae	NA	0.964
bin_Aruba.57	Actinomycetota	Actinomycetia	Mycobacteriales	Mycobacteriaceae	Corynebacterium pyruviciproducens	0.931
bin_Aruba.51	Actinomycetota	Coriobacteriia	Coriobacteriales	Atopobiaceae	Parolsenella sp900544655	0.969
bin_Denmark.44	Actinomycetota	Coriobacteriia	Coriobacteriales	Atopobiaceae	UBA7748 sp900314535	0.951
bin_Aruba.14	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	1.000
bin_Denmark.4	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	1.000
bin_Aruba.21	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	NA	0.999
bin_Brazil.61	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella tanakaei	0.998
bin_Aruba.39	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp902362275	0.992
bin_Spain.84	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp900555555	0.992
bin_Brazil.151	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella intestinalis	0.991
bin_Denmark.127	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp900548365	0.991
bin_Denmark.17	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella bouchesdurhonensis	0.983
bin_Brazil.81	Bacillota	Bacilli	RFN20	CAG-826	NA	0.907
bin_Brazil.109	Bacillota_A	Clostridia	Oscillospirales	Butyricicoccaceae	AM07-15 sp003477405	0.936
bin_Aruba.28	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.995
bin_Malaysia.9	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Dysosmobacter welbionis	0.986
bin_Malaysia.116	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Flavonifractor plautii	0.981
bin_Malaysia.34	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Lawsonibacter sp000177015	0.954
bin_Aruba.54	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.959
bin_CaboVerde.35	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.930
bin_Aruba.36	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	Peptoniphilus_C sp902363535	0.913
bin_Aruba.52	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.906
bin_Malaysia.26	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	NA	0.976
bin_Aruba.29	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	Gemmiger variabilis_C	0.941
bin_Malaysia.103	Bacillota_C	Negativicutes	Acidaminococcales	Acidaminococcaceae	Acidaminococcus intestini	0.942
bin_Malaysia.8	Bacillota_C	Negativicutes	Veillonellales	Dialisteraceae	Dialister sp900543165	0.981
bin_Denmark.60	Bacillota_C	Negativicutes	Veillonellales	Dialisteraceae	Allisonella histaminiformans	0.938
bin_Aruba.64	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera sp000417505	0.970
bin_Brazil.58	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera elsdenii	0.969
bin_CaboVerde.38	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera sp000417505	0.961
bin_Malaysia.81	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera stantonii	0.943
bin_Malaysia.96	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Caecibacter sp003467125	0.919
bin_Brazil.163	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella lascolaii	0.992
bin_CaboVerde.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella copri	0.981
bin_Aruba.10	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900544825	0.978
bin_Brazil.5	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900540415	0.929
bin_Spain.21	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola coprophilus	0.929
bin_Malaysia.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp000437675	0.924
bin_Malaysia.117	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp900541335	0.919
bin_Malaysia.77	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp900542985	0.913
bin_Brazil.75	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	NA	0.964
bin_Malaysia.61	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	Helicobacter_C labetoulli	0.964
bin_Brazil.128	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	NA	0.963
bin_Brazil.70	Desulfobacterota	Desulfovibrionia	Desulfovibrionales	Desulfovibrionaceae	Mailhella sp003150275	0.944
bin_Brazil.146	Pseudomonadota	Alphaproteobacteria	RF32	CAG-239	CAG-495 sp001917125	0.969
bin_Brazil.9	Pseudomonadota	Alphaproteobacteria	RF32	CAG-239	CAG-495 sp000436375	0.919
bin_Aruba.15	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.963
bin_Brazil.51	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	Sutterella wadsworthensis_A	0.963
bin_Spain.43	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	CAG-521 sp000437635	0.959
bin_Brazil.64	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.948
bin_Brazil.92	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	CAG-521 sp000437635	0.948
bin_Brazil.15	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.930
bin_CaboVerde.33	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum sp900543125	0.995
bin_Brazil.82	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum succiniciproducens	0.990
bin_CaboVerde.55	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum_A thomasii	0.939

8.3.2 Negatively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum,class,family)%>%
    tt()

genome	phylum	class	order	family	species	value
bin_Denmark.27	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus_B hirae	0.000
bin_Brazil.12	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus faecalis	0.001
bin_Brazil.170	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	NA	0.014
bin_CaboVerde.24	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus_E cecorum	0.020
bin_CaboVerde.16	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Limosilactobacillus reuteri	0.000
bin_Denmark.56	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Levilactobacillus brevis	0.000
bin_Malaysia.72	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Limosilactobacillus reuteri	0.000
bin_CaboVerde.60	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus agilis	0.001
bin_Malaysia.127	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus agilis	0.001
bin_Denmark.61	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Latilactobacillus sakei	0.003
bin_Malaysia.35	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus animalis	0.003
bin_Malaysia.4	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus ruminis	0.035
bin_CaboVerde.14	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus acidophilus	0.067
bin_Aruba.18	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus lutetiensis	0.001
bin_Brazil.22	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus pasteurianus	0.001
bin_Denmark.117	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Lactococcus garvieae	0.001
bin_Denmark.113	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus alactolyticus	0.002
bin_Brazil.69	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	Clostridium_A leptum	0.097
bin_Denmark.93	Bacillota_A	Clostridia	Lachnospirales	Anaerotignaceae	Anaerotignum sp001304995	0.074
bin_CaboVerde.3	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium_P perfringens	0.000
bin_Denmark.42	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium_P perfringens	0.000
bin_Brazil.136	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium thermobutyricum	0.041
bin_Brazil.3	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900555395	0.026
bin_Brazil.89	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900550235	0.034
bin_Denmark.63	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerostipes hadrus	0.040
bin_Denmark.19	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900556835	0.052
bin_Denmark.118	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Dorea_A longicatena	0.067
bin_Malaysia.108	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.081
bin_Brazil.113	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerobutyricum sp900754855	0.095
bin_Denmark.50	Bacillota_A	Clostridia	Peptostreptococcales	Peptostreptococcaceae	Peptacetobacter sp900539645	0.053
bin_Brazil.107	Bacillota_A	Clostridia	Monoglobales_A	UBA1381	CAG-41 sp900066215	0.016
bin_Brazil.159	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900541465	0.040
bin_Brazil.140	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_A sp900543175	0.046
bin_Malaysia.6	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900545035	0.088

8.3.3 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()

Rows: 202 Columns: 2
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): phylum, colors

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.4 Correlations

#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)

# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
        keep.tip(., tip=m$spNames) 


#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
          + (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean

matrix <- toPlot %>% 
      as.data.frame() %>%
      rownames_to_column(var="genome1") %>%
      pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
      mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
      mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
      ggplot(aes(x = genome1, y = genome2, fill = cor)) +
            geom_tile() + 
            scale_fill_gradient2(low = "#be3e2b",
                       mid = "#f4f4f4",
                       high = "#b2b530")+
            theme_void()

htree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

vtree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#create composite figure
grid.arrange(grobs = list(matrix,vtree),
             layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1)))

8.5 Predict responses

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("domestic","feral")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "origin", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(origin=rep(gradient,1000)) %>%
            pivot_longer(!origin,names_to = "genome", values_to = "value")

# weights:  9 (4 variable)
initial  value 101.072331 
final  value 91.392443 
converged

8.5.0.1 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

8.5.1 Positive associated functions at element level

# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0103	0.008804968	6.661144e-04	0.016709154	0.920	0.080	Biosynthesis	Nucleic acid biosynthesis	UDP/UTP
D0504	0.004538946	4.053932e-04	0.009749750	0.917	0.083	Degradation	Amino acid degradation	Methionine
D0507	0.003748020	2.758746e-05	0.007466919	0.901	0.099	Degradation	Amino acid degradation	Leucine
D0906	0.003766165	1.962339e-04	0.008167397	0.926	0.074	Degradation	Antibiotic degradation	Fosfomycin
D0208	0.010264085	2.155175e-03	0.018260856	0.942	0.058	Degradation	Polysaccharide degradation	Mixed-Linkage glucans
D0205	0.012642911	2.073360e-03	0.023229274	0.939	0.061	Degradation	Polysaccharide degradation	Pectin

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0219	-0.003901767	-0.008826791	-2.848873e-04	0.043	0.957	Biosynthesis	Amino acid biosynthesis	GABA
B0214	-0.021208689	-0.038785864	-3.609382e-03	0.072	0.928	Biosynthesis	Amino acid biosynthesis	Glutamate
B0204	-0.015449214	-0.031753343	-6.766342e-04	0.092	0.908	Biosynthesis	Amino acid biosynthesis	Serine
B0302	-0.004719002	-0.009783858	-6.819933e-04	0.034	0.966	Biosynthesis	Amino acid derivative biosynthesis	Betaine
B0310	-0.012814580	-0.023534635	-3.361749e-03	0.046	0.954	Biosynthesis	Amino acid derivative biosynthesis	Tryptamine
B0303	-0.011870796	-0.021137564	-2.779849e-03	0.060	0.940	Biosynthesis	Amino acid derivative biosynthesis	Ectoine
B0309	-0.007506723	-0.015190011	-1.212630e-04	0.098	0.902	Biosynthesis	Amino acid derivative biosynthesis	Putrescine
B0804	-0.016705505	-0.030010083	-3.381213e-03	0.053	0.947	Biosynthesis	Aromatic compound biosynthesis	Dipicolinate
B0603	-0.016546534	-0.031262013	-2.676829e-03	0.069	0.931	Biosynthesis	Organic anion biosynthesis	Citrate
B0601	-0.009185987	-0.017495258	-7.435595e-04	0.088	0.912	Biosynthesis	Organic anion biosynthesis	Succinate
B0401	-0.012323164	-0.023303496	-7.038749e-04	0.085	0.915	Biosynthesis	SCFA biosynthesis	Acetate
B0709	-0.002087118	-0.003665291	-6.444408e-04	0.032	0.968	Biosynthesis	Vitamin biosynthesis	Tocopherol/tocotorienol (E)
D0517	-0.004416113	-0.007674836	-1.257010e-03	0.026	0.974	Degradation	Amino acid degradation	Ornithine
D0508	-0.003149520	-0.007229581	-1.762986e-04	0.079	0.921	Degradation	Amino acid degradation	Lysine
D0903	-0.003862111	-0.008824495	-2.684042e-04	0.043	0.957	Degradation	Antibiotic degradation	Cephalosporin
D0908	-0.015343133	-0.027690198	-4.284755e-03	0.052	0.948	Degradation	Antibiotic degradation	Macrolide
D0601	-0.009191565	-0.017397410	-2.465709e-03	0.037	0.963	Degradation	Nitrogen compound degradation	Nitrate
D0610	-0.003027935	-0.004879338	-9.417130e-04	0.039	0.961	Degradation	Nitrogen compound degradation	Methylamine
D0611	-0.003862111	-0.008824495	-2.684042e-04	0.043	0.957	Degradation	Nitrogen compound degradation	Phenylethylamine
D0603	-0.001850442	-0.003639849	-3.827215e-04	0.048	0.952	Degradation	Nitrogen compound degradation	Urate
D0606	-0.005646982	-0.010629987	-1.125038e-03	0.060	0.940	Degradation	Nitrogen compound degradation	Allantoin
D0612	-0.001583928	-0.002951867	-1.645286e-05	0.098	0.902	Degradation	Nitrogen compound degradation	Hypotaurine
D0801	-0.001879207	-0.002299701	-1.080150e-04	0.004	0.996	Degradation	Xenobiotic degradation	Toluene
D0802	-0.001879207	-0.002299701	-1.080150e-04	0.004	0.996	Degradation	Xenobiotic degradation	Xylene
D0807	-0.004059376	-0.008574328	-6.102883e-04	0.055	0.945	Degradation	Xenobiotic degradation	Catechol
D0817	-0.004725492	-0.010236141	-6.085526e-04	0.056	0.944	Degradation	Xenobiotic degradation	Trans-cinnamate
D0816	-0.005551332	-0.011507872	-4.487159e-04	0.083	0.917	Degradation	Xenobiotic degradation	Phenylacetate

8.5.2 Negatively associated functions at element level

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0219	-0.003901767	-0.008826791	-2.848873e-04	0.043	0.957	Biosynthesis	Amino acid biosynthesis	GABA
B0214	-0.021208689	-0.038785864	-3.609382e-03	0.072	0.928	Biosynthesis	Amino acid biosynthesis	Glutamate
B0204	-0.015449214	-0.031753343	-6.766342e-04	0.092	0.908	Biosynthesis	Amino acid biosynthesis	Serine
B0302	-0.004719002	-0.009783858	-6.819933e-04	0.034	0.966	Biosynthesis	Amino acid derivative biosynthesis	Betaine
B0310	-0.012814580	-0.023534635	-3.361749e-03	0.046	0.954	Biosynthesis	Amino acid derivative biosynthesis	Tryptamine
B0303	-0.011870796	-0.021137564	-2.779849e-03	0.060	0.940	Biosynthesis	Amino acid derivative biosynthesis	Ectoine
B0309	-0.007506723	-0.015190011	-1.212630e-04	0.098	0.902	Biosynthesis	Amino acid derivative biosynthesis	Putrescine
B0804	-0.016705505	-0.030010083	-3.381213e-03	0.053	0.947	Biosynthesis	Aromatic compound biosynthesis	Dipicolinate
B0603	-0.016546534	-0.031262013	-2.676829e-03	0.069	0.931	Biosynthesis	Organic anion biosynthesis	Citrate
B0601	-0.009185987	-0.017495258	-7.435595e-04	0.088	0.912	Biosynthesis	Organic anion biosynthesis	Succinate
B0401	-0.012323164	-0.023303496	-7.038749e-04	0.085	0.915	Biosynthesis	SCFA biosynthesis	Acetate
B0709	-0.002087118	-0.003665291	-6.444408e-04	0.032	0.968	Biosynthesis	Vitamin biosynthesis	Tocopherol/tocotorienol (E)
D0517	-0.004416113	-0.007674836	-1.257010e-03	0.026	0.974	Degradation	Amino acid degradation	Ornithine
D0508	-0.003149520	-0.007229581	-1.762986e-04	0.079	0.921	Degradation	Amino acid degradation	Lysine
D0903	-0.003862111	-0.008824495	-2.684042e-04	0.043	0.957	Degradation	Antibiotic degradation	Cephalosporin
D0908	-0.015343133	-0.027690198	-4.284755e-03	0.052	0.948	Degradation	Antibiotic degradation	Macrolide
D0601	-0.009191565	-0.017397410	-2.465709e-03	0.037	0.963	Degradation	Nitrogen compound degradation	Nitrate
D0610	-0.003027935	-0.004879338	-9.417130e-04	0.039	0.961	Degradation	Nitrogen compound degradation	Methylamine
D0611	-0.003862111	-0.008824495	-2.684042e-04	0.043	0.957	Degradation	Nitrogen compound degradation	Phenylethylamine
D0603	-0.001850442	-0.003639849	-3.827215e-04	0.048	0.952	Degradation	Nitrogen compound degradation	Urate
D0606	-0.005646982	-0.010629987	-1.125038e-03	0.060	0.940	Degradation	Nitrogen compound degradation	Allantoin
D0612	-0.001583928	-0.002951867	-1.645286e-05	0.098	0.902	Degradation	Nitrogen compound degradation	Hypotaurine
D0801	-0.001879207	-0.002299701	-1.080150e-04	0.004	0.996	Degradation	Xenobiotic degradation	Toluene
D0802	-0.001879207	-0.002299701	-1.080150e-04	0.004	0.996	Degradation	Xenobiotic degradation	Xylene
D0807	-0.004059376	-0.008574328	-6.102883e-04	0.055	0.945	Degradation	Xenobiotic degradation	Catechol
D0817	-0.004725492	-0.010236141	-6.085526e-04	0.056	0.944	Degradation	Xenobiotic degradation	Trans-cinnamate
D0816	-0.005551332	-0.011507872	-4.487159e-04	0.083	0.917	Degradation	Xenobiotic degradation	Phenylacetate

positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

table <- bind_rows(positive,negative) %>%
  left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) 

table %>%
  mutate(Element=factor(Element,levels=c(table$Element))) %>%
  ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.5.2.1 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D02	8.380877e-03	-0.0032042156	0.021012695	0.829	0.171
B08	8.092933e-03	-0.0036562752	0.019064905	0.803	0.197
B01	1.198216e-03	-0.0057372614	0.008344889	0.652	0.348
S01	9.353796e-04	-0.0114173781	0.013508236	0.586	0.414
B09	4.002654e-05	-0.0005537729	0.000499940	0.361	0.639

# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D08	-0.0011456784	-0.002162438	-0.0002132115	0.039	0.961
B03	-0.0106204008	-0.018053351	-0.0033897581	0.050	0.950
D06	-0.0029947487	-0.006665266	0.0001147733	0.107	0.893
B04	-0.0084917327	-0.018498150	0.0007309851	0.130	0.870
B06	-0.0067968210	-0.016392244	0.0029246318	0.177	0.823
D07	-0.0121698848	-0.031193068	0.0045376311	0.192	0.808
D03	-0.0041778567	-0.013047116	0.0032739106	0.207	0.793
D05	-0.0014533798	-0.007604949	0.0034585469	0.224	0.776
S03	-0.0094794343	-0.032363817	0.0155818115	0.240	0.760
B02	-0.0034965044	-0.012889226	0.0053410419	0.278	0.722
D09	-0.0019870507	-0.008375102	0.0042295505	0.302	0.698
S02	-0.0046167761	-0.014546536	0.0029897655	0.332	0.668
B07	-0.0034403466	-0.016055313	0.0096265812	0.348	0.652
D01	-0.0001590696	-0.005078656	0.0040657195	0.434	0.566
B10	-0.0000088381	-0.000324572	0.0002401476	0.494	0.506

positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")